
#|________________________________________________________________________
 |
 | pcamain.lsp
 |
 | Principal Components ViSta Plugin Model Object
 | Copyright (c) 1991-2000 by Forrest W. Young 
 |
 | This file contains the plugin constructor function, 
 | defproto function and isnew method along with
 | slot accessor methods and system methods
 |
 | This code is loaded at first use time
 |_____________________________________________________________________
 |#


;;------------------------------------------------------------------------
;;define function for constructing a principal component plugin model object
;;------------------------------------------------------------------------

(defun principal-components 
  (&key 
     (data       *current-data*)
     (dialog      nil)
     (covariances nil))
"ViSta plugin function to perform Principal Components analysis.  With no arguments, calculates correlations among all active numeric variables in the current data and performs a principal components analysis of those variables. Keyword arguments are:
:COVARIANCES followed by 
   t    causes the analysis to be performed on covariances, or 
   nil  causes the analysis to be performed on correlations, the default.
:DATA   followed by the data-object to be analyzed (default: *current-data*)
:DIALOG followed by t (to display parameters dialog box) or nil (default)
[constructor function]"
  (send pca-model-object-proto :new *pca-menu-item-name* data dialog (not covariances)))

;;------------------------------------------------------------------------
;;define principal component model prototype object and its slot-accessors
;;------------------------------------------------------------------------

(defproto pca-model-object-proto 
  '(scores coefs eigenvalues svd corr corcovmat 
           var-rel-contrib obs-rel-contrib) 
  () vista-analysis-plugin-object-proto)

(defmeth pca-model-object-proto :isnew (*pca-menu-item-name* data dialog corr)
  (cond
    ((> (send current-data :active-nvar '(numeric)) 
        (send current-data :active-nobs))
     (error-message (format nil "Note: Cannot analyze data ~
                with fewer active observations (~d) ~
                than active numeric variables (~d)." 
                        (send current-data :active-nobs) 
                        (send current-data :active-nvar '(numeric))))
     (send *toolbox* :reset-button tool-name))
    ((< (send current-data :active-nvar '(numeric)) 2)
     (error-message (format nil "Cannot analyze data ~
                that has less than two active numeric variables."))
     (send *toolbox* :reset-button tool-name))
    (t
     (send self :corr corr)
     (call-next-method *pca-menu-item-name* data dialog))))


(defmeth pca-model-object-proto :scores (&optional (values nil set))
  (if set (setf (slot-value 'scores) values))
  (slot-value 'scores))

(defmeth pca-model-object-proto :coefs (&optional (values nil set))
    (if set (setf (slot-value 'coefs) values))
    (slot-value 'coefs))

(defmeth pca-model-object-proto :eigenvalues (&optional (values nil set))
    (if set (setf (slot-value 'eigenvalues) values))
    (slot-value 'eigenvalues))

(defmeth pca-model-object-proto :svd (&optional (structure nil set))
  (if set (setf (slot-value 'svd) structure))
  (slot-value 'svd))

(defmeth pca-model-object-proto :corr (&optional (val nil set))
  (if set (setf (slot-value 'corr) val))
  (slot-value 'corr))

(defmeth pca-model-object-proto :corcovmat (&optional (matrix nil set))
  (if set (setf (slot-value 'corcovmat) matrix))
  (slot-value 'corcovmat))

(defmeth pca-model-object-proto :var-rel-contrib (&optional (matrix nil set))
  (if set (setf (slot-value 'var-rel-contrib) matrix))
  (slot-value 'var-rel-contrib))

(defmeth pca-model-object-proto :obs-rel-contrib (&optional (matrix nil set))
  (if set (setf (slot-value 'obs-rel-contrib) matrix))
  (slot-value 'obs-rel-contrib))

(defmeth pca-model-object-proto :options ()
  (when (send self :dialog)
        (let ((result nil)
              (dialog-value (choose-item-dialog
                "Analysis Options:"
              '("Analyze Correlations"
                "Analyze Covariances") :initial 0))
              )
          (when dialog-value
                (when (= 0 dialog-value) (setf result t))
                (send self :corr result))
          dialog-value)))

(defmeth pca-model-object-proto :analysis ()
  (let* (;(left-alpha 1)
         (data (send self :data-matrix))
         (nobsnvar (array-dimensions data))
         (nobs (select nobsnvar 0))
         (nvar (select nobsnvar 1))
       #|  (prepped-data 
           (if (send self :corr) 
               (/ (normalize (center data) 1) (sqrt (1- nobs)))
               (/ (center data) (sqrt (1- nobs)))))|#
         (prepped-data   (if (send self :corr) 
                             (normalize (center data))
                             (center data)))
         (
         (corcovmat (covariance-matrix prepped-data))
         (svd (sv-decomp2 prepped-data)) 
         (svd (if (< (sum (col (select svd 2) 0)) 0)
                     (list (* -1 (select svd 0))
                           (select svd 1)
                           (* -1 (select svd 2))
                           (select svd 3))
                  svd))
         (scores (matmult (select  svd 0) 
                          (diagonal (^ (select svd 1) left-alpha))))
         (eigenvalues (^ (select svd 1) 2))
         (coefs (matmult (select  svd 2) 
                         (diagonal (^ (select svd 1)(1- left-alpha)))))
         (var-rel-contrib (^ (matmult coefs (diagonal (select svd 1))) 2))
         (sqrt-obs-coord-row-list (row-list (^ (* scores (sqrt nobs)) 2)))
         (obs-rel-contrib (apply (function bind-rows) 
                                 (/ sqrt-obs-coord-row-list 
                                    (mapcar #'sum sqrt-obs-coord-row-list))))
         )
    (send self :nobs nobs)
    (send self :nvar nvar)
    (send self :corcovmat corcovmat)
    (send self :svd svd)
    (send self :coefs coefs)
    (send self :scores scores)
    (send self :eigenvalues eigenvalues)
    (send self :var-rel-contrib var-rel-contrib)
    (send self :obs-rel-contrib obs-rel-contrib)
    t))


(defmeth pca-model-object-proto :fitted (i ndim)
"Arg: I NDIM
Computes fitted data values for variable I given an NDIM dimensional fit (e.g., using the first NDIM componnents scores and coefficients."
  (let* ((nobs (send self :nobs))
         (nvar (send self :nvar))
         (scores (select (send self :scores) (iseq nobs) (iseq ndim))) ;nobs x ndim
         (evals (diagonal (select (send self :eigenvalues) (iseq ndim)))) ;ndim x ndim
         (coefs (select (transpose (send self :coefs)) (iseq ndim) (list i)))) ;ndim x 1
    (matmult scores (^ evals 2) coefs)))


(defmeth pca-model-object-proto :save-model-template (data-object)
"Args: (data-object)
DATA-OBJECT is the object-identification information of a data object. 
The method contains a template for saving the model-object." 
  `(principal-components    
    :title      ,(send self :title)
    :name       ,(send self :name) 
    :dialog      nil
    :covariances ,(not (send self :corr))
    :data (data  ,(send data-object :name)
                 :title      ,(send data-object :title)
                 :variables ',(send self :variables)
                 :types     ',(send self :types)
                 :labels    ',(send self :labels)
                 :data      ',(send self :data)))
  )

(defmeth pca-model-object-proto :create-data 
  (&key (dialog nil)
        (scores t)
        (coefs  t set-coefs)
        (coords nil)
        (rel-contribs nil)
        (abs-contribs nil)
        (input  nil))
"Args: DIALOG (SCORES T) COEFS COORDS REL-CONTRIBS ABS-CONTRIBS INPUT
Creates output data objects. If DIALOG=T then presents dialog to determine which objects created. Otherwise presents specified objects. If no options, specified, creates only SCORES."
  (if (not (eq current-object self)) (setcm self)) 
  (let ((creator (send *desktop* :selected-icon))
        (desires (list (list (if scores 0) 
                             (if coefs 1)
                             (if coords 2)
                             (if rel-contribs 3)
                             (if abs-contribs 4)
                             (if input 5)
                             )))
        )
    (cond 
      (dialog
       (setf desires 
             (choose-subset-dialog "Choose Desired Data Objects"
                  '("Component Scores"
                    "Component Coefficients"
                    "Coordinates"
                    "Relative Contributions"
                    "Absolute Contributions" 
                    "Analyzed Input Data")
                   :initial (select desires 0))))
      (t
       (setf desires 
             (list (list (if scores 0)
                         (if (and set-coefs  coefs ) 1)
                         (if coords 2)
                         (if rel-contribs 3)
                         (if abs-contribs 4)
                         (if input 5))))))
    
    (when desires
          (when (member '0 (select desires 0))
                (send current-model :pca-scores-data-object creator))
          (when (member '1 (select desires 0))
                (send current-model :pca-coefs-data-object  creator))
          (when (member '2 (select desires 0))
                (send current-model :pca-coords-data-object  creator))
          (when (member '3 (select desires 0))
                (send current-model :pca-rel-contrib-data-object  creator))
          (when (member '4 (select desires 0))
                (send current-model :pca-abs-contrib-data-object  creator))
          (when (member '5 (select desires 0))
                (send current-model :create-input-data-object "PCA" creator)))
    (not (not desires))))

(defmeth pca-model-object-proto :pca-scores-data-object (creator)
  (data (strcat "Scores-" (send self :name))
   :created creator
   :creator-object self
   :title (strcat "PCA Scores for " (send self :title))
   :data (combine (send self :scores))
   :variables (mapcar #'(lambda (x) (format nil "PC~a" x)) 
                      (1+ (iseq (min (send self :nvar) (send self :nobs)))))
   :labels (send self :labels)
   :types (repeat "Numeric" (min (send self :nvar) (send self :nobs)))
  ))

(defmeth pca-model-object-proto :pca-coefs-data-object (creator)
  (data (concatenate 'string "Coefs-" (send self :name))
   :created creator
   :creator-object self
   :title (concatenate 'string "PCA Coefficients (Eigenvectors) for " 
                       (send self :title))
   :data (combine (send self :coefs))
   :labels (send self :variables)
   :variables (mapcar #'(lambda (x) (format nil "PC~a" x)) 
                      (1+ (iseq (send self :nvar))))
   :types (send self :types))
  )

(defmeth pca-model-object-proto :pca-coords-data-object (creator)
  (let* ((coefs (send self :coefs))
         (eigenvalues (send self :eigenvalues))
         (var-coords (matmult coefs (diagonal (sqrt eigenvalues))))
         (scores (send self :scores))
         (nobs (send self :nobs))
         (obs-coords (* scores (sqrt nobs)))
         (nvar (send self :nvar))
         )
    (data (concatenate 'string "Var-Coords-" (send self :name))
          :created creator
          :creator-object self
          :title (concatenate 'string "PCA Variable Coordinates for " 
                              (send self :title))
          :data (combine var-coords)
          :labels (send self :variables)
          :variables (mapcar #'(lambda (x) (format nil "PC~a" x)) (1+ (iseq nvar)))
          :types (send self :types))
    (data (concatenate 'string "Obs-Coords-" (send self :name))
          :created creator
          :creator-object self
          :title (concatenate 'string "PCA Observation Coordinates for " 
                              (send self :title))
          :data (combine obs-coords)
          :variables (mapcar #'(lambda (x) (format nil "PC~a" x)) 
                             (1+ (iseq (min nvar nobs))))
          :labels (send self :labels)
          :types (repeat "Numeric" (min nvar nobs)))
    ))

(defmeth pca-model-object-proto :pca-rel-contrib-data-object (creator)
  (let* ((var-rel-contrib (send self :var-rel-contrib))
         (obs-rel-contrib (send self :obs-rel-contrib))
         (nobs (send self :nobs))
         (nvar (send self :nvar))
         )
    (data (strcat "Var-Rel-Contrib-" (send self :name))
          :created creator
          :creator-object self
          :title (strcat "PCA Relative Contributions of Variables for " 
                         (send self :title))
          :data (combine var-rel-contrib)
          :labels (send self :variables)
          :variables (mapcar #'(lambda (x) (format nil "PC~a" x)) (1+ (iseq nvar)))
          :types (send self :types))
    (data (strcat "Obs-Rel-Contrib-" (send self :name))
          :created creator
          :creator-object self
          :title (strcat "PCA Relative Contributions of Observations for " 
                              (send self :title))
          :data (combine obs-rel-contrib)
          :variables (mapcar #'(lambda (x) (format nil "PC~a" x)) 
                             (1+ (iseq (min nvar nobs))))
          :labels (send self :labels)
          :types (repeat "Numeric" (min nvar nobs)))
    ))

(defmeth pca-model-object-proto :pca-abs-contrib-data-object (creator)
  (let* ((coefs (send self :coefs))
         (eigenvalues (send self :eigenvalues))
         (var-coords (matmult coefs (diagonal (sqrt eigenvalues))))
         (scores (send self :scores))
         (var-abs-contrib 
          (matmult (^ var-coords 2) (diagonal (/ 1 eigenvalues))))
         (obs-abs-contrib
          (matmult (^ scores 2) (diagonal (/ 1 eigenvalues))))
         (nobs (send self :nobs))
         (nvar (send self :nvar))
         )
    (data (strcat "Var-Abs-Contrib" (send self :name))
          :created creator
          :creator-object self
          :title (strcat "PCA Absolute Contributions of Variables for " 
                         (send self :title))
          :data (combine var-abs-contrib)
          :labels (send self :variables)
          :variables (mapcar #'(lambda (x) (format nil "PC~a" x))  
                             (1+ (iseq nvar)))
          :types (send self :types))
    (data (strcat "Obs-Abs-Contrib-" (send self :name))
          :created creator
          :creator-object self
          :title (strcat "PCA Absolute Contributions of Observations for " 
                         (send self :title))
          :data (combine obs-abs-contrib)
          :variables (mapcar #'(lambda (x) (format nil "PC~a" x)) 
                             (1+ (iseq (min nvar nobs))))
          :labels (send self :labels)
          :types (repeat "Numeric" (min nvar nobs)))
    ))


(defmeth pca-model-object-proto :report (&key (dialog nil))
  (if (not (eq current-object self)) (setcm self))   
  (let* ((corcov (send self :corr))
         (corcovst (if corcov "Correlations" "Covariances"))
         (choices 
          (if dialog
              (choose-subset-dialog "PRNCMP Report Options:"
                                    (list corcovst
                                          "Coordinates"
                                          "Contributions"
                                          ))
              (list nil)))
         (w nil) (w1 nil) (w2 nil) (w3 nil) (yesno t)
         (corcovst nil) (labels nil) (vars   nil)
         (scores nil) (coefs  nil) (eigenvalues nil)
         (proportions nil) (coordinates nil) (fitmat nil) (lc-names nil)
         (nobs (send current-data :active-nobs))
         )
    (when choices
          (setf w (report-header (send self :title) :page t))
          (setf corcovst (if corcov "Correlation" "Covariance"))
          (setf labels (send self :labels))
          (setf vars   (send self :variables))
          (setf scores (send self :scores))
          (setf coefs  (transpose (send self :coefs)))
          (setf eigenvalues (send self :eigenvalues))
          (setf proportions (/ eigenvalues (sum eigenvalues)))
          (setf coordinates (matmult (transpose coefs) 
                                     (diagonal (sqrt eigenvalues))))
          (setf fitmat (transpose 
                        (matrix (list 3 (min (send self :nobs) 
                                             (send self :nvar))) 
                                (combine eigenvalues proportions
                                         (cumsum proportions)))))
          (setf lc-names (mapcar #'(lambda (x) (format nil "PC~a" x)) 
                                 (1+ (iseq (send current-model :nvar)))) )
          (when (= (length choices) 1)
                (when (not (first choices)) (setf choices nil)))
          (display-string 
           (format nil "Principal Components Analysis of Variable") w)
          (display-string (format nil " ~a ~2%"corcovst) w)
          (display-string (format nil "Model:     ~a~%" (send self :name )) w)
          (display-string (format nil "Variables: ~a~%" vars) w)
          (when choices
                (when (member 0 (first choices))
                      (display-string (format nil "~%~a Matrix~%"corcovst) w)
                      (print-matrix-to-window (send self :corcovmat) w 
                                              :row-heading "Variables"
                                              :column-heading "Variables"
                                              :row-labels vars
                                              :col-labels vars
                                              :decimals 4)))
          (display-string (format nil "~%") w)
          (display-string (format nil "Fit Measures for each Component:~%Eigenvalue (amount of total data variance fit by each component)~%Proportion (of total data variance fit by each component)~%Cumulative Proportion (of total data variance fit by the components)~2%") w)
          (print-matrix-to-window fitmat w 
                                  :row-heading "Components"
                                  :column-heading "Fit Measures"
                                  :col-labels '("E-Value" "Prop." "CumProp")
                                  :row-labels lc-names 
                                  :decimals 5)
          (display-string (format nil "~%Coefficients (Eigenvectors):~%") w)
          (print-matrix-to-window (transpose coefs) w 
                                  :row-heading "Variables"
                                  :column-heading "Components"
                                  :col-labels lc-names
                                  :row-labels vars
                                  :decimals 4)
          (display-string (format nil "~%Component Scores:~%(Left Singular Vectors times Square Root of Eigenvalues)~%") w)
          (print-matrix-to-window scores w
                                  :row-heading "Observations"
                                  :column-heading "Components"
                                  :col-labels lc-names
                                  :row-labels labels :decimals 4)
          (send w :scroll 0 0)
          (when choices
                (when (member 1 (first choices))
                      
                      (display-string (format nil "~4%COORDINATES:~2%Variable Coordinates:~%(Eigenvectors times Square-Root of Eigenvalues)~%") w)
                      (print-matrix-to-window 
                       coordinates w 
                       :row-heading "Variables"
                       :column-heading "Component Coordinates"
                       :row-labels vars 
                       :column-labels lc-names
                       :decimals 4)
                      (display-string (format nil "~%Observation Coordinates:~%(Scores times Square-Root of number of observations)~%") w)
                      (print-matrix-to-window 
                       (* scores (sqrt nobs)) w 
                       :row-heading "Observations"
                       :column-heading "Component Coordinates"
                       :column-labels lc-names
                       :row-labels labels 
                       :decimals 4)
                      (send w :scroll 0 0))
                (when (and (member 2 (first choices))
                           (if (> (send self :nobs) 200)
                               (yes-or-no-dialog 
                                (format nil "NOTE: The output for contributions may be long and~%slow. Do you really want to see this?"))
                               t))
                    
                      
                      (display-string 
                       (format nil "~4%CONTRIBUTIONS:~2%Absolute Contribution of Variables:~%(Squared Coordinates divided by Eigenvalues)~%") w)
                      (print-matrix-to-window
                       (matmult (^ coordinates 2) 
                                (diagonal (/ 1 eigenvalues))) 
                       w 
                       :row-heading "Variables"
                             :column-heading "Component Contributions"
                       :column-labels lc-names
                             :row-labels vars 
                       :decimals 4)
                      (display-string
                       (format nil"~%Absolute Contribution of Observations:~%(Squared Scores divided by Eigenvalues)~%") w)
                      (print-matrix-to-window 
                       (matmult (^ scores 2) (diagonal (/ 1 eigenvalues)))
                       w 
                       :row-heading "Observations"
                       :column-heading "Component Contributions"
                       :column-labels lc-names
                       :row-labels labels 
                       :decimals 4)
                      (display-string
                       (format nil "~%Relative Contribution of Variables:~%(Squared Coordinates)~%") w)
                      (print-matrix-to-window 
                             (^ coordinates 2) w 
                       :row-heading "Variables"
                       :column-heading "Component Contributions"
                       :column-labels lc-names
                       :row-labels vars 
                       :decimals 4)
                      (display-string
                       (format nil "~%Relative Contribution of Observations:~%(Squared Coordinates divided by the row-sum of squares)~%") w)
                      (let* ((sqrt-obs-coord-row-list 
                              (row-list (^ (* scores (sqrt nobs)) 2)))
                             (what-you-want
                              (/ sqrt-obs-coord-row-list
                                 (mapcar #'sum sqrt-obs-coord-row-list)))
                             (matrix-of-what-you-want
                              (apply (function bind-rows) what-you-want))
                             )
                        (print-matrix-to-window 
                         matrix-of-what-you-want w 
                         :row-heading "Observations"
                         :column-heading "Component Contributions"
                         :column-labels lc-names
                         :row-labels labels 
                         :decimals 4))
                      (send w :scroll 0 0))
          ))
  (send w :fit-window-to-text)
  w))    
  
(defun scree-plot (eigenvalues &key location size)
"Args: (eigenvalues)
Takes an eigenvalue sequence, produces a scree-plot, returns a plot-object."
   (let* ((n (length eigenvalues))
          (prop (/ eigenvalues (sum eigenvalues)))
          (scree-plot (plot-points (+ 1 (iseq n)) prop
                :show nil :title "Scree Plot"
                :variable-labels '("Principal Components" "Proportion")))
          (maxy (select prop 0))
          (cum-prop (cumsum prop)))
     (send scree-plot :add-lines 
           (list (+ 1 (iseq n)) prop) :draw nil :color 'red)
     (send scree-plot :point-color (iseq n) 'red)
     (mapcar #'(lambda (i)
                 (send scree-plot :point-label i
                       (cond
                         ((> i 0) 
                          (format nil "~5,4f, ~5,4f, ~5,4f" 
                                  (select prop i) 
                                  (- (select prop (- i 1)) (select prop i) )
                                  (select cum-prop i)))
                         (t
                          (format nil "~5,4f" (select prop i))))))
             (iseq (send scree-plot :num-points)))
     (send scree-plot :range 1 0  (* .1 (ceiling (* 10 maxy))) :draw nil)
     (send scree-plot :y-axis t t (1+   (ceiling (* 10 maxy))) :draw t)
     (send scree-plot :showing-labels t)
     (send scree-plot :mouse-mode 'brushing)
     (send scree-plot :menu nil)
     (send scree-plot :point-state 0 'selected)
     (send scree-plot :adjust-to-data :draw nil)

     (defmeth scree-plot :plot-help ()
       (plot-help-window (strcat "Help for " (send self :title)))
       (paste-plot-help (format nil "The Scree plot shows the relative fit of each principal component. It does this by plotting the proportion of the data's variance that is fit by each component versus the component's number. The plot shows the relative importance of each component in fitting the data.~2%"))
       (paste-plot-help (format nil "The numbers beside the points provide information about the fit of each component. The first number is the proportion of the data's variance that is accounted for by the component. The second number is the difference in variance from the previous component. The third number is the total proportion of variance accounted for by the component and the preceeding components.~2%"))
       (paste-plot-help (format nil "The Scree plot can be used to aid in the decision about how many components are useful. You use it to make this decision by looking for an elbow (bend) in the curve. If there is one (and there often isn't) then the components following the bend account for relatively little additional variance, and can perhaps be ignored.~2%"))
       (show-plot-help))
     scree-plot))

(require (strcat *prncmp-plugin-constructor-file* "pcavis"))
(require (strcat *prncmp-plugin-constructor-file* "vision"))
(require (strcat *prncmp-plugin-constructor-file* "vision+"))